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Abstract 

Theories of solvation free energies often involve electrostatic potentials at the 
position of a solute charge. Simulation calculations that apply cutoffs and 
periodic boundary conditions based on molecular centers result in center- 
dependent contributions to electrostatic energies due to a systematic sorting 
of charges in radial shells. This sorting of charges induces a surface-charge 
density at the cutoff sphere or simulation-box boundary that depends on the 
choice of molecular centers. We identify a simple solution that gives correct, 
center-independent results, namely the radial integration of charge densities. 
Our conclusions are illustrated for a Lennard-Jones solute in water. The 
present results can affect the parameterization of force fields. 
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Accurate simulation calculations of free energies of solvation require a careful treat- 
ment of long-range electrostatic interactions. Recent computational and theoretical work 
on single-ion free energies^ has converged upon a common set of ideas that are, however, 
discussed in slightly different ways, i.e., Gaussian fluctuations of electrostatic potentials^ 



the calculation of electrostatic potentials at atom positions on a solute molecule at frac- 
tional charge states {e.g., uncharged or fully charged). However, a lack of consensus on how 
electrostatic potentials should be evaluated means that calculated partial contributions to 
single-ion free energies are often not fully comparable. Differences arise because of a common 
practice of evaluating electrostatic interactions considering whole molecules. This can lead 
to spurious dependences on the choice of the center of a molecule. Similar issues arose in 
calculations of the electrostatic potential difference of the water-vapor interface: seemingly 
identical calculations of electrostatic potentials can produce different final results.! 

Discrepancies in calculated electrostatic potentials were noted recently by Aqvist and 
Hansson.i The present paper resolves the difficulties noted there. We will focus on the 
calculation of electrostatic potentials at the position of a solute molecule in a polar fluid, 
discussing the effects of different methods of summing charge interactions. This leads us to a 
simple, center independent, and feasible recipe used to analyze electrostatic potentials, both 
in finite and infinite systems, namely spherical integration of charge densities. To illustrate 
our general results, we will show data for Lennard- Jones (LJ) solutes in water. 

Two different center dependences will be considered (see Figure [I]). The first is associated 
with the center of the solvent molecule denoted by M used to bin electrostatic interactions be- 
tween solvent molecules and the solute molecule. The second center dependence to consider 
is the dependence on the solvent center P that might be used in implementing minimum- 
image periodic boundary conditions (PBC's) by translating a whole solvent molecule into 
the primary simulation box. These two centers M and P would often coincide but they need 
not. The effects considered are distinct. 

For molecule-based summation, the electrostatic potential at the center of a spherical 



second-order perturbation theory,@ or linear-response 




approaches require 
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LJ solute molecule depends strongly on the choice of the center M of a water molecule 
that defines into which shell it belongs. Shown in Figure are curves 0m ( r ) of potential 
contributions of water molecules with their center M within a radius r of the solute molecule, 

Mr) = r dr /£] 5(r - r hM ) J2 —) , (1) 

where the a sum extends over the water oxygen atom O and hydrogen atoms HI and H2. 
(. . .) denotes a canonical ensemble average over a system of iV SPC water molecules^ with 
oxygen and hydrogen positions r i; o, Yi,m and r^H2, respectively {r iy o = \?i,o\, etc.); and one 
uncharged LJ solute atom at position r s = with SPC-water LJ parameters. 5(r) is the 
Dirac delta function. q Q and qu are the charges on the oxygen and hydrogen sites of SPC 
water (— 0.82e and 0.41e, respectively). r it M is the center of water molecule i, defined as 
r i,M = wr it o + (1 — w)(r iim + r i)H2 )/2. The atom positions r it o, r iim and r ijH2 are shifted 
molecule-based under PBC's. (That is, the center P=M is mapped into the simulation box, 
leaving the molecule intact so that individual atoms can actually be outside the simulation 
box.) For weights w — 1 and 0, the center position r ijM coincides with the oxygen position 
and the hydrogen bisector, respectively. 

The molecule-based potential defined in eq |l| contrasts with the charge-based potential 
<f) q (r): 



4> q (r) = Arc r dr p q (r)/r , (2a) 



o 
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PjT) = EE ^(47rr 2 )- 1 5(r - r* J . (2b) 



\i=l a=l 



p q (r) is the radially averaged charge density. In eq ^T|, PBC's for the positions of charges 
r| G , rf Hl and rf H2 are applied on the basis of atoms rather than molecules. 

Each of the <pAi(r) curves in Figure |2] for different centers M reaches a plateau value 
after 0.6 to 0.8 nm distance from the solute. However, the plateau values differ not only in 
magnitude but also in sign for different choices of M, whereas identical choices of M give 
agreement between simulations under PBC's and using clusters with 256 and 1024 water 
molecules. The differences are caused by the M-dependent sorting of molecules, even for 



identical configurations (positions and orientations) of the solvent molecules. If the center 
M is close to the oxygen atom, the first layer of molecules considered in the integration 
in eq [I] predominantly includes water molecules with the oxygen atoms facing the solute. 
Correspondingly, 4>m{x) starts out negative as negative contributions of the oxygen atoms 
dominate. On the other hand, if the center M is close to the hydrogen atoms, the first 
layer of molecules considered in the integration will predominantly have the hydrogen atoms 
facing the solute (see Figure [j], middle). As a consequence, (puij) starts out positive and 
also reaches a positive plateau value. Results for centers M between the oxygen and the 
hydrogen bisector fall between the two curves. 

For a finite sample, the different curves all converge to the same value when all contribu- 
tions have been summed up (Figure |2j, middle and bottom). Convergence is therefore reached 
only after crossing the interface to the exterior, so that surface-potential contributions are 
included. For the cluster simulations of Figure (middle and bottom), the potential crosses 
a liquid- vacuum interface. 

Similar problems arise with molecule-based cutoffs (or residue-based cutoff's for macro- 
molecules). For instance, if the distance to the oxygen atom of a water molecule is used to 
determine whether a particle interacts with that water molecule, a characteristic surface- 
charge density is induced at the cutoff sphere. The oxygen density seen by the solute is 
essentially a step function. The hydrogen density is reduced just inside the cutoff and 
nonzero just beyond the cutoff, resulting in a net negative charge density just inside the 
cutoff sphere and a net positive charge density just outside. This effective surface-dipole 
density strongly affects the potential at the site of the particle. That effect is independent 
of the cutoff length, as the surface area and charge-dipole interaction vary with the square 
and the inverse square of the cutoff length, respectively. 

When whole molecules are shifted under PBC's this leads to another level of ill-definition 
of electrostatic potentials. Shifting molecules as a whole means that PBC's are applied based 
on a center P with coordinates r i)P = ur i + (1 — u)(r i H1 + r itH2 )/2 with weight u. If that 
center P coincides with the center M of the <pM integration, then the plateau in ^Af(r) 
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reached after the first layer of water molecules remains essentially unchanged when reaching 
the box boundary. However, if P and M do not coincide, 0m (r) crosses over from the M 
curve to the P curve when the box boundary is reached. This can be seen in Figure [2] (top) 
where the <f>M(f) curve for M equal to O (w — 1, u — 0) crosses over to the hydrogen-bisector 
curve (w = 0,u = 0) when the hydrogen-bisector is used for PBC's (P=HH). Clearly, this 
is an unphysical behavior associated with summing electrostatic interactions and applying 
PBC's on the basis of molecules. 

How can we eliminate these difficulties of calculating electrostatic potentials in computer 
simulations? The unphysical false plateaus observed for <pM {r) in Figure |] stem from asso- 
ciating partial charges with molecular centers. By choosing a center M, the water molecules 
were systematically sorted for analysis. For a finite system, integration to infinity is re- 
quired to get the correct result. And that result will then contain troublesome and undesired 
surface-potential contributions. Under PBC's, that integration cannot be performed easily, 
as is manifest from the dependence of the limiting value of the potential on the choice of 
the molecular center P upon which PBC's are applied. 

However, if we alternatively integrate over charge densities p q (v) rather than sum over 
molecules, we will obtain a well-defined result for the potential that coincides with taking 
the limit of an infinite system before extending the integral to infinity. The charge-based 
potential <j> q (r) is defined in eq |a|. For a finite system, eqs [l] and || will give identical results 



if the integration volume covers the whole system (extending beyond the interface to the 
container, vacuum etc.). However, unlike eq |]the potential 4> q (r) defined in eq|2a| will reach 
a plateau beyond the correlation length of the charge correlation p q (r) independent of an 
arbitrary choice of the center M of a molecule. (As shown in Figure |^, that plateau is reached 
within about 1 nm from the neutral LJ solute. Larger correlation lengths were observed for 
a charged solute.i) 

These issues would be largely irrelevant with conventional Ewald treatment of electro- 
static potentials, where the simulation box is replicated periodically in space. However, 
center dependences can arise with modifications of the standard Ewald approach. The 



electrostatic potentials of periodic images can be summed up using the Ewald potential 
V?s(r).i0 ^Pe(t) is the periodic solution of Poisson's equation A<£#(r) = — Air[5(r) — l/V] 
for a unit point charge and a homogeneous background in the unit cell V. The equivalents 
of the electrostatic potentials 0m (f) and (j> q (r) defined in eqs [L| and |2a] for periodic systems 
are then 

r I N 3 V 

0mO) = / dr(J2 6 ( r ~ t ^m) H Qa<fiE(r i>a ) j , (3a) 

J ° \i=l a=l I 

r I N 3 \ 

tff(r) = I dr ( £ £ % - ^ ) • (3b) 

Again, minimum-image PBC's for and r| a are applied on the basis of molecular centers 
P and individual atoms, respectively. Figure [3| shows that the charge-based Ewald potential 
and 1/r curves <ftq(r) and <p q {r) converge but that the molecule-based curve 4>M{r) for 
periodic systems also converges to 4>q{r) rather than (puij"). This is expected because the 
Ewald potential is fully periodic. 

Physical modification of the Ewald potential sacrifice this periodicity. The Ewald poten- 
tial is the limit of performing the lattice sum with the growing lattice embedded in a sphere 
cut out of a medium with infinite dielectric constant e' = oo (tin-foil boundary conditions). 
Total potential energies without the effect of that dielectric background e' = oo require sub- 
traction of a term proportional to the square of the net dipole moment M of the simulation 
box.0 Expressed as an effective potential, we can subtract a term 27rr 2 /3l / from ^(r): 
<fE,e'=i( r ) — ¥e{ y ) — 2nr 2 /3V. This destroys the periodicity. Use of the modified potential 
<fE,e'=i( r ) m ec l |H forces 0f/(r) to converge to 0Af( r ), as shown in Figure |3|. However, the 
result for the potential (pfjir) at the solute site then again depends on the particular choice 
of the molecular center P upon which PBC's are applied. Clearly, to reproduce the non- 
physical effects of integrating the potential using 1/r with molecule-based sorting requires 
subtraction of a non-periodic term from the Ewald potential and application of the potential 
outside the "universe," i.e., the simulation box. 

It must be noted that subtracting the r 2 term from the Ewald potential has little effect 
if the integration is based on charges (Figure^). However, applying the r 2 modification by 
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molecules rather than atoms leads to large differences in the potential. The charge-based 
potential (eq |3ap with P=0 molecular-based shifting of the r 2 term has a value different in 
sign from the charge-based curves shown in Figure We also emphasize that changing 
the dielectric background to a finite value e' < oo in the Ewald sum should not affect the 
charging of an ion at the center of the box. The dipolar field induced by a background 
e' beyond a spherical cavity around r = is proportional to r • M which is zero at the 
position r = of the uncharged particle. When a point multipole is charged from zero, that 
contribution is also zero because of averaging over all orientations. 

The results of this paper explain the differences in the sign of the electrostatic potential 
at the position of an uncharged LJ particle in water between Aqvist and Hanssoni (M=0 
based sorting, 1/r: negative potential), Rick and Berne0 (charge-based sorting; Ewald 
and r 2 modification with P=0 based shifting: negative potential) and Pratt at alM as 
well as Hummer et alM (charge based sorting; Ewald, 1/r and a generalized reaction-field 
interaction: positive potential). The best current value for that potential is positive. In 
that context a re-examination of several results regarding free energies of charged species 
might be worthwhile. For instance, free energies of anions were found to be less negative in 
Ref . ^ than in Ref . [5] but more negative for cations. That can be explained if molecule-based 



summation has been used in Ref. [L3| using a center M at or close to the oxygen atom of 
water. The present results also affect the parameterization of force fields involving charged 
species. Finally, we emphasize that the errors induced by molecule-based summation are 
independent of the cutoff length for sufficiently large cutoffs. If the induced surface- charge 
distribution were symmetrically distributed on a spherical shell then it follows from Gauss's 
law that the correction to the induced electrostatic potential inside the spherical shell would 
be a constant. In that case, the contributions of M-dependent sorting would cancel each 
other for an overall neutral, polar solute but not for a solute with a net charge. 

Our results suggest that these issues are primarily matters of analysis of configurational 
simulation data. A variety of methods may be used to obtain the configurational data. The 
center dependences considered here are introduced by the analysis of electrostatic poten- 



tials and are often larger than the secondary differences in the configurational data due to 
variations in their production. 

The following general recipe for electrostatic-potential calculations emerges: (1) Elec- 
trostatic interactions should be integrated based on charge densities rather than individual 
molecules to give correct results for atoms and molecules carrying point charges or spatially 
extended charge distributions. For molecule-based summation, the calculated potentials 
4>(r) level out nicely but the plateau values depend on the arbitrary choice of molecular 
centers. (2) In simulations using PBC's, all charges should be mapped into the simulation 
box. Molecule-based PBC's result in center-dependent surface-charge densities. (3) Un- 
der PBC's, Ewald summation provides an accurate way of summing up all interactions, 
minimizing finite-size effects. 

Note added in proof. Ashbaugh and WoodS come to similar conclusions regard- 
ing molecule-center dependences of electrostatic potentials in their comparison of Ewald 
summation! and cutoff calculations.^ In particular, these authors also find the potential to 
be positive for a neutral LJ solute in water. 
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FIGURES 

FIG. 1. M- and P-center sorting of molecular partial charges. Left: Different M-centers con- 
sidered for the water molecule. M and M' coincide with the hydrogen bisector and the oxygen 
position. Middle: Charges of the upper molecule are counted in the shaded spherical shell (bin) 
but not the charges of the lower molecule. The lower molecule with an outward-pointing dipole 
moment is placed in a more distant bin. Right: P-center sorting, where P coincides with the 
oxygen position. The bottom image of the molecule is considered in the electrostatic potential 
calculations. For the particular choice of P=0 and isotropic molecular orientations, the charge 
density is depleted just inside the simulation cell around the solute (outlined as square and circle, 
respectively) and enriched just outside. 

FIG. 2. Integrated electrostatic potentials at the position of an uncharged LJ solute in SPC 
water using 1/r interactions. Results are shown for different ways of sorting the charges and 
applying PBC's (atom or molecule based). The top panel shows the results of averaging over 
140 000 Monte-Carlo passes of a system with 255 SPC water molecules and one LJ solute with 
SPC-water LJ parameters (using Ewald summation; see Ref. ^ for simulation details). M and P 
denote the centers of sorting and applying PBC's, respectively, where O is the oxygen and HH the 
hydrogen-bisector position. The middle and bottom panel show the results of averaging over 100 000 
and 300 000 Monte-Carlo passes of clusters of 256 and 1024 SPC water molecules, respectively, and 
one LJ particle at the center, again with SPC-water LJ parameters. In the cluster simulations, 
electrostatic interactions were calculated using 1/r Coulomb interactions without cutoff. The 
asymptotic value of charge-based integration using the Ewald potential is shown for reference. 

FIG. 3. Integrated electrostatic potential at the position of an uncharged LJ solute in SPC water 
using the Ewald potential <pe( v ) instead of 1/r. Results are shown for charge and molecule-based 
integration with and without the r 2 modification added to <pe(v). See Figure |2| for further details. 
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